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Using highly efficient simulations of the tight-binding Bogoliubov-de Gennes model we solved 
self-consistently for the pair correlation and the Josephson current in a Superconducting-Bilayer 
graphene-Superconducting Josephson junction. Different doping levels for the non-superconducting 
link are considered in the short and long junction regime. Self-consistent results for the pair cor- 
relation and superconducting current resemble those reported previously for single layer graphene 
except in the Dirac point where remarkable differences in the proximity effect are found as well as a 
suppression of the superconducting current in long junction regime. Inversion symmetry is broken 
by considering a potential difference between the layers and we found that the supercurrent can be 
switched if junction length is larger than the Fermi length. 

PACS numbers: 73.43.-f, 73.23.-b, 73.63.-b 



I. INTRODUCTION 

Ever since graphene, a two dimensional crystal with 
honeycomb lattice structure, became recently available 
in research laboratories, a new variety of hybrid struc- 
tures has been explored. This lead to the development 
of promising technological devices and the understanding 
of the non-intuitive physical mechanisms of rclativistic- 
like massless fermions. It is well-known that electrons 
propagating through single layer graphene (SLG) exhibit 
many peculiar properties, e.g. a gapless linear dispersion 
around the neutrality (or Dirac) point which resemble 
massless Dirac fermions, a minimum conductivity at zero 
carrier concentration, high mobility as well as the absence 
of backscattering. When in contact to superconductors, 
graphene exhibits exotic superconducting properties. Al- 
though graphene was not found to sustain intrinsic su- 
perconductivity itself, there is experimental evidence 1 - - — 
that when in proximity with a conventional superconduc- 
tor it becomes superconducting. Because of the conven- 
tional superconducting proximity effect, which describes 
how Cooper pairs diffuse from the superconducting ma- 
terial into metals, superconducting hybrid structures like 
superconducting-normal-superconducting (SNS) Joseph- 
son junctions rise as an interesting systems in which 
one can study electronic correlation of relativistic-like 
particles. These experiments have attracted consider- 
able theoretical attention insofar as it was predicted that 
a finite supercurrent should exist even at zero doping 
where the density of states is vanishing^ or that specular 
Andreev reflections should happen at the superconduc- 
tor/graphene interface in SLG 7 . In the ballistic regime, 
which is realizable in SNS graphene Josephson junctions 
according to recent experiments^, theoretical studies 
have shown the existence of a finite bipolar supercon- 
ducting current through the junction^ - — . 

With a gapless parabolic, instead of linear, band struc- 
ture, bilayer graphene (BLG) appears at this point as a 
suitable alternative for investigating electronic correla- 
tions in two dimensional systems. Also, as it is widely 
known theoretically and experimentally, a tunable gap 



can be induced in BLG by an out-of-plane applied electric 
field, which is very useful for transistor application s 13 ' 14 . 
We expect that the role of the gapless parabolic disper- 
sion may be important since the superconducting corre- 
lations depend strongly on the electronic properties of 
the material. Moreover, the ability to open a gap in the 
spectrum by an external electric field could be of interest 
in superconducting devices. 

Within a tight-binding Bogoliubov-de Gennes formal- 
ism we calculated self-consistently the pair correlation 
and the Josephson current through a Superconducting- 
BLG-Superconducting Josephson junction. Our findings 
resemble previous numerical results shown for SLG-based 
Josephson junctions except near the Dirac point in case 
of a long junction where a suppression of the current 
is found, mainly due to the vanishing density of states 
present only in non-dimer sites. We also show how the 
superconducting current can be switched off by applying 
a out-of-plane electric field. 



II. MODEL 

We consider bilayer graphene in the common AB 
(Bernal) stacking with two inequivalent sublattices A and 
B in the top layer, and the corresponding sublattices B 
and A in the bottom layer. We model the interlayer inter- 
action by a hopping parameter which couples the nearest- 
neighbors in sub-lattices A and B from the top and bot- 
tom layer. Any additional interlayer coupling terms are 
ignored in this study. It is well known that in this case the 
electron dispersion displays two parabolic bands touching 
at the Dirac point and two additional parabolic bands at 
higher energies due to interlayer induced splittin g 13 : 15 . 

We study the Josephson effect in a BLG-based junc- 
tion by considering the hybrid nanostructure illustrated 
in Fig. [1] The top and bottom layers are both in con- 
tact with superconducting leads while a junction of size 
L is suspended. The influence of the superconducting 
contacts is modeled by assuming an on-site attractive 
pairing potential, Us < and heavy doping, /is > in 
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layers, makes it necessary to perform the self-consistent 
calculation for all four sub-lattice types within the unit 
cell. 

In order to describe the DC Josephson effect we fix the 
difference in the phases of the order parameter, A<fi s — 
4>R — <t>L, in the outermost parts of the S regions (see 
Fig. [T]). In a finite region near the interfaces, with size 
of the order of the coherence length, the order parameter 
is allowed to relax self-consistently. Therefore the phase 
gradient over the non-superconducting region A(j) N will 
be restricted according to Acj) N < A(j) S < tt, considering 
that the maximum value of A(j) S is tt. This constraint for 
A^ was recently pointed out by Black-Shaffer et al~2. 
in a SLG-based ballistic junction. 



FIG. 1. Layout of the SNS-BLG Josephson junction where 
the superconducting leads are modeled by assuming on-site 
attractive pairing potential Us and a heavy doping [is in 
the regions under the contacts labeled by S. For the non- 
superconducting link, with length L and labeled by TV, we 
took the pairing potential Un = and a varying chemical 
potential /ijv. Phases 4>l and 4>R are kept fixed during the 
self-consistently calculation. 



the contact regions, which are labeled by S, such that 
a s-wave superconducting state is induced in the outer- 
most regions of both layers. The normal region labeled 
by TV, which has a tunable Fermi level [In and zero pair- 
ing potential, Un — 0, acts like a non-superconducting 
channel through which Cooper pairs could tunnel. Simi- 
lar models for graphene based Josephson junctions were 
previously considered in the Dirac limit&i^r— as well as 
in the tight-binding formulatio n 12 i 16 i 17 . We solve self- 
consistently for the order parameter along the junction 
following closely the self-consistent calculation performed 
in Ref. [l2| for a ballistic single-layer graphene Josephson 
junction. This is necessary in order to consider the possi- 
bility of Cooper pairs being depleted close to the interface 
in the superconducting region due to the existence of the 
normal region, i.e. the inverse proximity effect. As is usu- 
ally assumed for SLG-based junctions we consider a clean 
and smooth interface such that all physical quantities are 
homogeneous along the i-direction parallel to the inter- 
face. In addition, we have considered a wide junction, 
with width W » L, as well as periodic boundary condi- 
tions imposed along the y-direction. The latter assump- 
tions make it possible to reduce the self-consistent cal- 
culation of the order parameter from a two-dimensional 
problem to a one-dimensional one since we can restrict 
the calculation to only a unit cell along the y-direction 
perpendicular to the interface. It is worth mentioning 
that the pairing potential does not distinguish between 
sublattices defined on either layer. However, due to the 
interlayer coupling, the self-consistent calculation should 
be performed separately for only two inequivalent sites 
A and B (or A and B) within the unit cell. The break- 
ing of the inversion symmetry, which is easily achieved 
by considering a potential difference applied between the 



III. NUMERICAL METHOD: CHEBYSHEV 
EXPANSION OF THE GREEN'S FUNCTION 

Superconducting correlations in a BLG-based Joseph- 
son junction are described by using the following tight- 
binding Hamiltonian with on-site attractive Hubbard in- 
teractions: 



U 



( 



- Yl + ei ) 4*^ + + e2 ) 4^ w 

i 

where ctjvac) creates a spin-up electron on the i-site in 
the top layer whereas Cjjjvac) creates a spin-down hole on 
the j-site in the bottom layer. The hopping parameter, 
t, describes the intralayer hopping integral between next- 
nearest neighbors in the same layers while t± = 0.143t 
correspond to the interlayer nearest-neighbors hopping 
which couples the dimer sites A and B. Other hopping 
terms, like the interlayer coupling between the non-dimer 
B and A sites, are not considered in the present work, 
since they influence only very low-energy excitations. 
The Fermi level is shifted from the charge neutrality point 
or Dirac point by the chemical potential \ii and Ui is the 
on-site attractive pairing potential which is non vanish- 
ing only in the right and left superconducting regions. 
The on-site energies t\ and ti for atomic sites on the top 
and bottom layer, respectively, have been introduced in 
order to simulate a potential difference or gate voltage 
V g = t\ — £2 between the layers. 

By using the Hartree-Fock decomposition and keep- 
ing only terms relevant to the superconducting order, 
one can transform the many-body Hamiltonian |T]) into a 
mean-field single-particle Hamiltonian, which within the 
Nambu formalism can be written as follows: 
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where Ho and A are the following 2x2 matrices: 



considering that all physically quantities calculated here, 
such as density of states, pair correlation or the Joseph- 
son current, can be obtained from the Green's function 
for each lattice point separately. 
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A, 
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where the diagonal elements of the matrix Q correspond 
to the on-site mean-field superconducting order parame- 
ter A l = Ui(cif ai). 

Following Refs. LL8| and LL9| we have performed the 
self-consistent mean-field calculation through a numeri- 
cal approximation of the Gorkov Green's function by us- 
ing the Chebyshev-Bogoliubov-de-Gennes method. Both, 
the normal and anomalous Gorkov Green's function, can 
be approximated by a superposition of a finite number of 
Chebyshev polynomials as follows: 
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where the expansion coefficients for the diagonal, or nor- 
mal (a — 1), and the off-diagonal, or anomalous (a = 2), 
components of the 2x2 Green function are defined re- 
spectively a$*£: 
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where T n (x) — acos(n cos(x)) is the Chebyshev polyno- 
mial of order n, which is defined according to the recur- 
rence relation: T n+ i(x) — 2xT{x) — T n -\(x). 

In order to perform the sum ([5]) one needs first to 
rescalc the Hamiltonian such that the eigenvalues lie 
in the [-1,1] interval. To this end, matrices ([2]) and (HJ as 
well as the energies have being normalized according to: 
H -> H = {H- lb) /a andw^w = (w - lb) /a, where 
the rescaling factors are a = (E max — E m i n )/{2 — rj) and 
b = (E max + E m i n )/2, with rj > being a small number. 

Once the Hamiltonian is normalized, the expansion co- 
efficients can be obtained through a recursive procedure. 
Starting with an initial vector | jo) = |c^) and the first 
term in the iteration = | jo) one can obtain the 
n-th term by using the Chebyshev recurrence relation: 
\j n ) = 2H\j n -i) — \jn-2}- Chebyshev moments are fi- 
nally obtained from the scalar product (a\j n ), where (a\ 
are the vectors (1| = {c^\ and (2| = (c|, | for moments 
expanding the (O and ([7]) component of the Green's func- 
tion respectively. 

Since most of the computational effort corresponds 
to sparse matrix-vector and vector-vector multiplication, 
high speed-up can be achieved by implementing parallel 
computation on graphical processing units (GPUs). We 
are therefore able to solve efficiently systems described 
by matrices of sizes between 88Kx88K and 320Kx320K 
according to the junction size considered in this study. 
Additional parallel computations can be implement by 



IV. RESULTS 

It is well-known that the energy dispersion in SLG and 
BLG differ around the neutrality point. Therefore, qual- 
itative differences are expected in the proximity effect as 
well as in the Josephson current in the two systems. For a 
quantitative comparison between SLG and BLG Joseph- 
son junction we have set up the following values for the 
physical input parameters: Us — — 1.36t = — 3.4eV and 
fis = 0.6t = 1.5eV, similar to the values used in a pre- 
vious self-consistent study for ballistic SLG Josephson 
junctions^. These parameters lead to a finite s-wave 
bulk superconducting order parameter, Ao = 0.041i, 
which corresponds to a superconducting coherence length 
£ = hvp/Af) s» 23 unit cells similar to the one consid- 
ered in the SLG junction case. Due to the difference 
in the local density of states between dimer and non- 
dimer locations, the order parameter is slightly differ- 
ent for these two types of atoms. Both junction length 
regimes are solved with the proposed self-consistent nu- 
merical method, e.g. short junction for L < £ and long 
junction for L > £. 

Previous analytical descriptions of BLG Josephson 
junctions based on the Dirac equatio n 10 ' 20 requires 
smooth interfaces and a low energy regime, for which 
Ao << t± << /is. Here these restrictions are lifted but 
in order to compare with relevant experimental scenarios 
we performed calculations only for situations correspond- 
ing to A < tj_ < MS- 



A. Proximity effect 

We show in Fig. [5] the self-consistently calculated pair 
correlation (c^c^) for both inequivalent sites, A (dimer) 
and B (non-dimer), in a unit cell defined along the in- 
direction perpendicular to the interface. As was pre- 
viously mentioned, both long and short junctions are 
considered here and plotted in Figs. [5Ja) and[^b), re- 
spectively. We present the pair correlation function for 
several doping levels in the non-superconducting region 
showing that the proximity effect is strongly dependent 
on the relative Fermi level mismatch between the S and 
N regions. A significant difference between BLG and 
SLG can be seen for the undoped case (jUjv = 0) while 
for other dopings the pair correlation is found to exhibit 
similar behavior in BLG and SLG. In particular, we can 
see that in the undoped case dimer sites in BLG show 
a suppression of the pair correlation over the N region 
compared to the SLG case. Opposite behavior is seen 
for non-dimer sites where a larger pair leakage into N is 
found. This behavior is similar to the proximity effect 
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FIG. 2. Absolute value of the pair correlation in a long (a) and 
short (b) BLG Josephson junction as a function of the position 
along the y-direction perpendicular to the SN interface. Both 
inequivalent sites, dimer at A and non-dimer at B, are plotted 
separately for different doping levels considered for the non- 
superconducting region: y,ff — 0, 0.3t and 0.6t. The first one 
correspond to the case when the Fermi level is pinned at the 
Dirac point while the last one corresponds to no-FLM at the 
interface, (c) Pair correlation at L/2 in the N region as a 
function of the junction length L for different values of fiN- 
SLG self-consistent results are shown for comparison with the 
BLG undoped case. 



in strained graphene where a sublattice polarization of 
the local density of states in the zero-th pseudo-Landau 
level induces sub-lattice dependent leaking distances^. 
No relevant differences in the pair correlation profile are 
found between dimer and non-dimer sites in BLG for 
higher doping levels considered here: moderately doped 
/ijv = 0.3t and highly doped /ijv = 0.6i. Note also in Fig. 
[2Kb) that the reported intcrlaycr asymmetry in the pair 
correlation in BLG is found to be not important for the 
short junction regime. In fact, we can clearly see in Fig. 
[UJc) that the difference between the pair correlation at 
sites A and B becomes larger as the junction length is 
increased, as large as a few orders of magnitude. 
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FIG. 3. LDOS for A (red line) and B (blue dashed line) 
sub-lattice sites in the N region with (a) zero doping (at the 
Dirac point) and (b) heavy doping (or with no-FLM between 
the superconducting and normal regions), (c) LDOS in the 
highly doped S region showing the superconducting gap. 



superconducting regions. The LDOS is plotted for dimer 
(A) and non-dimer (B) lattice sites along the inhomoge- 
neous direction in the BLG Josephson junction away from 
the SN interface. Three particular cases have been cho- 
sen and shown in Fig. [3] Panel (a) shows the LDOS for 
two sites in the N region for the undoped case, whereas, 
in panel (b) we have depicted LDOS for the same loca- 
tion but when there is no-FLM. Results shown in Fig. 
|3Ja)-(b) are consistent with the fact that the density of 
states at dimer sites vanishes linearly around the Dirac 
point, while being finite at the non-dimer site o 13 ' 15 . 

Due to the differences in the LDOS for inequivalent 
sites in BLG the proximity effect will give different leak- 
ing distances in different sub-lattices, as we show in Fig. 
[5Jc) for the undoped case. Increasing the doping level in 
the normal region leads to the LDOS becoming homo- 
geneous. This result is consistent with the fact that no 
difference between both kind of sites is observed for the 
other doped cases in the pair correlation function shown 
in the last section. Finally, panel (c) in Fig. [3] shows the 
LDOS deep inside the superconducting region where we 
can see clearly the coherence peaks on each side of the su- 
perconducting gap. Note that in the normal region, the 
LDOS is modified near the Fermi level, corresponding to 
the formation of Andreev levels. 



C. Josephson current 



B. LDOS 

To further understand this peculiar behavior observed 
only in the undoped case we have plotted the local den- 
sity of states (LDOS) in the superconducting and non- 



1. Unbiased case 

In the absence of applied bias voltage but in the pres- 
ence of a finite phase difference between the two super- 
conducting sides, a DC supercurrent will flow across the 
junctio n 12 i 21 . This is the usual DC Josephson effect. For 
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FIG. 4. Self-consistent Josephson current as a function of 
position along the y-direction for top (a) and bottom (b) layer 
as well as the average current (c) in the system for different 
values of the pairing potential, i.e. different coherence lengths. 
The interlayer current between the A-B dimer sites is plotted 
in (d). 




FIG. 5. Current-phase relation normalized to the critical cur- 
rent found in (e) for SLG and BLG-based junctions consid- 
ering different doping levels in the N region: (e),(f) undoped 
/xjv=0, (c),(d) slightly doped /ijv=0.1t and (a),(b) no-FLM 
fJ-N=fJ-s=0.6t for short (a),(c),(e) and long (b),(d),(f) junc- 
tions. The short and long junction lengths are L = 10 and 
L — 50, respectively. 



this purpose a phase bias is achieved by fixing a desired 
phase difference between the outermost parts of the su- 
perconducting regions, <f>L and 4>r for the left and right 
sides of the junction respectively (see Fig. [T]). In or- 
der to numerically calculate the Josephson current we 
solve self-consistently for both phase and amplitude of 
the order parameter along the junction except in the 
extreme regions where we keep the phases fixed. The 
current profile along the junctions in both layers as well 
as the interlayer current are shown in Figs. IUa)-(d) f° r 
an undoped non-superconducting link for different values 
of the pairing potential. As we can see, the supercur- 
rent is found not to be constant within each layer sepa- 
rately, contrary to what is usually expected to happen 
for self-consistent current calculation in 2-dimensional 
systems- 1 . Instead, one interesting feature in BLG is 
the appearance of a weak interlayer current between A-B 
dimer sites as a consequence of the current conservation 
law. We observe that the LDOS at the left interface is 
asymmetric in top-bottom layers while for the right in- 
terface the asymmetry is reversed. Because of this, the 
current is enhanced at the left interface in the top layer 
while being suppressed in the bottom layer, therefore a 
weak interlayer current appears. The reverse happens at 
the right interface, but the average current remains flat 
across the whole junction, as expected. Next we con- 
struct the current-phase relation (CPR) by performing 
self-consistent current calculations for different phase dif- 
ferences between the superconducting contacts and dop- 
ing levels in the non-superconducting region. We show 
in Fig. [5] the phase-dependence of the current for both 
SLG and BLG for no-FLM cases with /in = 0.6t, panels 

(a) and (b), slightly doped with fi]y = 0. It < t± , panels 

(b) and (c) and undoped, panels (e) and (f) considering 



short, panels (a)-(c)-(e), and long junction, panels (b)- 
(d)-(f), regimes. Note that, in most of the cases shown 
in Fig. [5] a complete description of the current over the 
full [0,7r] phase-range is not possible. This restriction ap- 
pears as a consequence of the relaxation of the phase over 
the S region which becomes significant as the FLM is re- 
duced in the junction until S and N regions are equally 
doped and the phase drop goes linearly through the self- 
consistent region. Therefore the phase difference over 
the normal region will always be smaller than the ap- 
plied phase difference in the superconducting regions. A 
similar constraint is found for SLG-based junctions, as 
is shown in Fig. [5] and which was previously pointed out 
by Black-Schaffer el alH. We find that for no-FLM sit- 
uations the Josephson current density and the CPR of 
the BLG is almost identical to the one in SLG for both 
short and long junction. In contrast, for the undoped and 
slightly doped cases the short and long junctions have dif- 
ferent behavior. Short BLG junctions are similar to short 
SLG ones, while for long BLG junctions the Josephson 
current is suppressed. The origins of this suppression 
could be traced back to the sublattice polarization of 
the leaking distance which remains even for doping lev- 
els lower than the interlayer hopping energy, t±. While 
one sublattice (non-dimer sites) has an enhanced leak- 
ing distance, the other sublattice (dimer sites) behaves 
like an insulator with a short leaking distance. The re- 
sulting combination corresponds to a slightly suppressed 
Josephson current when compared to SLG. 
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FIG. 6. Critical density current in a short BLG junction, 
where L = 40 Unit cells considering a lower pairing potential 
Us = — 1.2t, as a function of the gate voltage Vg at the Dirac 
point (a) and the chemical potential /ijv (b). In the left figure 
two different regimes are considered: L < Xf and L > Xf, 
showing an enhancement and suppression of the supercon- 
ducting current, respectively. On right side /i-dependence is 
plotted for bias and unbias cases. SLG is included for com- 
parison. 



2. Biased case 

Inversion symmetry can be broken in the BLG nanos- 
tructure by considering a potential difference, V g — £1 — 62 
with ei = Vg/2 and £2 = — V g /2, applied between the 
layers. As a consequence, a tunable gap Ao is induced 
at the Dirac point for the undoped case and therefore 
inversion symmetry breaking appears as a good possibil- 
ity to switch off the superconducting current when the 
voltage induced gap overtakes the superconducting gap, 
A g > Ao. At this point, we only consider short junc- 
tions where such proposal might have potential techno- 
logical applications. In Fig. Eta) we show the depen- 
dence of the supercurrent on the applied bias voltage, 
V g , perpendicular to the BLG layer. We observe two 
distinct regimes, according to the relation between the 
Fermi wave-length and the junction length. In the case 
when £ > \ F = hvF/iVg/2) > L, we observe an enhance- 
ment of the Josephson current. This is because the charge 
density in the junction does not recover its bulk expected 
value when the BLG is under bias (+n for one layer and 
— n in the other layer) but has a finite positive value. In 
this regime the junction is effectively doped, thus showing 
an enhancement of the current. A different dependence 



of the current is achieved when £ > L > Xp. In this case 
for the center of the junction the charge density has the 
opposite polarity, and the expected suppression of the 
current with bias voltage is obtained. 

In Fig. Etb) we show the /^-dependence of the current 
for the biased and unbiased cases. A doping activation 
is found in the bias case for an energy value around A g , 
in addition a slight increase in the current is observed at 
jitjv ~ A g due to the enhancement of the LDOS above the 
gap edge. SLG results are included in order to compare 
directly with BLG and show that for the chosen length 
the current density is slightly larger in SLG. In the in- 
set of Fig. |U[b) we focus on the low doping regime and 
observe that in BLG additional oscillations appear, rem- 
iniscent of what was previously found analytically in the 
Dirac approximation^. 



V. CONCLUSIONS 

In conclusion by using an efficient numerical proce- 
dure we solve self-consistently the Bogoliubov de Gennes 
equations for a tight-binding model of the AB-stacked 
bilayer graphene Josephson junction. When compared 
to single layer graphene Josephson junctions we uncover 
several regimes. First, in the short junction regime, the 
current density is similar for SLG and BLG for any dop- 
ing of the normal junction region. Second, in the long 
junction regime, for undoped junctions, the BLG cur- 
rent density is suppressed while for doped junctions (with 
doping larger than t±) the BLG and SLG junction be- 
have in a similar way. We attribute the peculiar behav- 
ior of the undoped BLG junctions to the difference of the 
LDOS between the dimer and non-dimer sites, which give 
suppressed or enhanced Cooper pair leaking distances 
depending on the sublattice. We have calculated the 
current-phase relation and showed that similar to SLG, 
even for short junctions there is a departure from con- 
ventional symmetric CPR. Finally, we have shown that 
by applying a gate voltage perpendicular to the BLG a 
gap in the spectrum can be induced and a supercurrent 
switch can be achieved given that the junction length is 
larger than the Fermi wave-length. 
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